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Abstract 

The method of geodesic deviations has been applied to derive accurate analytic approximations to 
geodesies in Schwarzschild space-time. The results are used to construct analytic expressions for the 
source terms in the Regge- Wheeler and Zerilli-Moncrief equations, which describe the propagation of 
gravitational waves emitted by a compact massive object moving in the Schwarzschild background 
space-time. The wave equations are solved numerically to provide the asymptotic form of the wave 
at large distances for a series of non-circular bound orbits with pcriastron distances up to the ISCO 
radius, and the power emitted in gravitational waves by the extreme-mass ratio binary system is 
computed. The results compare well with those of purely numerical approaches. 
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1 Introduction 



In a recent paper [T] we have implemented a new kind of perturbation theory for geodesies 
in curved space-time, an improved version of the method of geodesic deviations [21 [3l HJ |5] . 
In contrast to the well-estabhshed post-newtonian scheme [6l[7], this perturbation theory is 
fully relativistic by construction, and is expected to become especially relevant for motion in 
strongly curved space-time regions. This holds in particular for matter moving in the vicinity 
of black-hole horizons. Indeed, our method seems particularly well-suited to describe extreme 
mass-ratio binary systems, formed by a compact object -e.g., a neutron star or stellar- mass 
black hole- orbiting a giant black hole such as found in the center of many galaxies, with 
masses exceeding a million solar masses and sometimes much more. 

Like the well-known compact binary systems of neutron stars [U [9], these extreme mass- 
ratio binaries are expected to lose energy by the emission of gravitational radiation. Discovery 
and analysis of this radiation would be a direct way to probe the strong-curvature region 
around the giant black hole. In this paper we summarize our covariant perturbation scheme 
and apply it to the calculation of gravitational waves emitted by extreme-mass ratio binaries. 
We will do this for the rotationally invariant case of Schwarzschild geometry, describing a non- 
rotating black hole. Gravitational waves are represented by the linear metric perturbations 
propagating on this background geometry, in particular the perturbations created by a test 
mass orbiting the black hole. Such perturbations are described by the Regge- Wheeler and 
Zerilli-Moncrief equations [TT] [T^ , the solutions of which represent the two physical modes 
of gravitational waves in a Schwarzschild background. 

Our perturbative solution for the geodesies allows us to write the source terms of the 
Regge- Wheeler-Zerilli-Moncrief equations in analytical form. The equations can then be 
solved numerically using the algorithm of Lousto and Price \13\ [T^ [T5] . Good agreement 
with the results existing in the literature is obtained. By just varying the initial data many 
solutions can be found in a very efficient way. 



2 Parametrization of orbits 



Complete exact solutions for the geodesic equations in Schwarzschild space-time are known 
for special cases, including circular orbits, the straight plunge and special inspiraling motions 

[T7] . More general orbits can be constructed in a covariant scheme as perturbations of 
these special analytic solutions; in this paper we consider especially bound orbits as covariant 
deformations of circular ones, using the results obtained in [1]. 

The orbits of test masses are described by geodesies x^^{t) parametrized by the proper time 
r. We choose standard Schwarzschild-Droste co-ordinates to parametrize the line element in 
the form 

\ _ 9M\ - - r'^de'' - sin2 . 



dr' = (^1 - — J - 2M - - ^ ^ d^ ■ (1) 

Moreover, as the conservation of angular momentum guarantees all orbits to be planar, we 
can fix the plane of the orbit to be the equatorial plane 6 = 7r/2. With initial conditions 
t(0) = 9^(0) = 0, circular orbits are then given by 

t{r)= r = R, ^(r) = \l ^ ^^=, (2) 

R 
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where R is the constant radial co-ordinate of the motion. Circular geodesies arc character- 
ized by special values of the constants of motion (e, £) corresponding to energy and angular 
momentum per unit of test mass: 



V r J dT A _ 3M dr Y 1 - 

For more general orbits e and i are unrelated, and r is not constant. Using the same initial 
conditions bound orbits can be parametrized as deformations of circular orbits of the form 

oo 

i(r) = OgT -|- sin nwr, 

n=l 

oo 

r(r) = + ^o; cos nwr, (4) 

n=l 
oo 

^{t) = OqT + aj^ sin ncor, 

n=l 

where both the amplitudes an and the angular frequency oj are to be expanded as power series 
in terms of a deformation parameter a, the series for the co-efficients an starting at order a". 
The parameter a is defined in terms of the distance between the perturbed and the circular 
orbit at r = 0, 

da'^ = g^ij dx^dx''\^^Q . (5) 

Analytic approximations to the orbits are obtained by cutting the expansions at some finite 
order in a: 



3M 
R 



flg = i? -I- crAi -I- ^ C7^A2 + a[ = an[ + ^a'^n2 + a2 = ^a'^ml + 

ip / 1 , f , ^ 2 , f f , 2 f , <fi ^ 2 f , 

^o=\ m I +avl + -a + a\ = anX + -a + a^ = -a + 

V -ft /i — 5M ^ II 

(6) 

with cocfRcicnts a{^, n > 3, and the dots representing terms of order cr'^ and higher. Also the 
angular frequency is to be computed order by order in the expansion parameter: 

w = cjo + ^'^i + ••• (7) 

The explicit expressions for the various terms above are given by 
a. for the time co-ordinate: 



t A 3M 1 

(A 2 9 36M I IS.SM'' 1&2M^ r2 i , M ' 

(l-^)(l-^) Rl-^-f 



(8) 
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and 

n\ = -2n{ xl ^ " = -Ani, 



6M 
R 



^ (5-^)(l-^)' 

r2 o 15M I 14M^ 
n ^ ^ ^ i ^2- 

2i? (1-^) 



4 = -A 



mn = A 



'2 

b. for the radial co-ordinate: 



2M 

IT' 



A r j\ O 15M _L 30M2 



-2 ^ (5-lf^)(l-^)(l-^)' 



7^2 



r2 1 _ 7M 
r _ "-1 R 



R 1-^ 



with A„, n = 1, 2, free parameters; 
c. for the angular co-ordinate: 



3Ai /m 1-^ 

A / R = Ai 

2i?2 V (1 _ 3M^3/2 - '^^1' 

/ a2i 25M I 153M^ 360M3 , 324M^ r2 i i M' 

I ^2 -I- 2M\ /i 3M\ {-i 6M^ p i 2M 



(11) 



2R (1-1^)(1-1^)(1-1^) 1- 



and 



^ 2nl I 

nl = — = -urii, 

R 

2niA, 1-^ + ir 

^2 — M 



R (5-^)(l-^)' 

r 2 _ 32M 
R . 
IM ' 



(12) 



4ii; 1 



"^2 — ^ n -, 6M 



d. and finally, for the angular frequency: 



I Ml- 'J 3Ai 1-1^ + 1^ 

""^-V^r^' - """^ ^ (1 - (1-^)- ^ ^ 

For fixed mass M, the undetermined constants in these expressions are (i?, c, Ai, A2, ...)• 
In every order there is one additional parameter to be fixed by initial conditions: R to fix 
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the circular reference orbit, a to fix the initial distance from the circular orbit (here: the 
periastron), and the shape parameters A„ determining the apastron, periastron advance, and 
other orbital characteristics. 

By construction the orbits Q imply the existence of the constants of motion e and i 
defined in eq. ([3|: 

\ r J dr' dr^ 

such that 

^=0, ^=0. (14) 

dr dr 

In terms of the expressions Q, ^ these identities indeed hold order by order in a for all 
allowed values of {R, a, A„), the restrictions being R > 6M and r > 2M. The values of e and 
£ are then given to order o"^ by 

e = £o + aei + ^a'^e2 + ..., i = £0 + (^h + ^(r'^h + (15) 



with 



_ 1 - ^ _ 1 oj^RAi 

^0 — / > ^1 — 7; I = 

'1 _ 3M 2 L 



R \ ^ R 

^ I ulRA, Mjnlf l-gf + ^r 

288A/ I 1279Af2 1836A/3 , 972Af4 



i?4 



and 





(17) 



243M I 963Af^ _ 1512Af'^ 972Af"t 
R ^ R^ R3 i;4 

Given values of the parameters (i?, cr, Ai, A„) up to some order n, we obtain a curve 
approximating a geodesic up to and including terms of order a"', and the corresponding 
values of the constants of motion e and i. Going to order n + 1 whilst keeping the values 
of the parameters (i?, a, Ai, .., A„) fixed then changes the approximate values of e and i, 
in a way depending on the choice of A„+i. In the limit n — t- 00 the curve approaches an 
exact geodesic, but in general this geodesic is not characterized by the initial values of the 
constants of motion for the circular orbit. The large redundancy introduced by the infinite 
set of parameters A„, which together determine the radius Oq of the fundamental circular 
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orbit, is easily understood, as any circular orbit can be deformed continuously into any other 
periodic (bound) orbit, including other circular orbits. The redundancy can be lifted in a 
simple way by taking all An = 0, which reduces the expansion to the one constructed in ref. 



[3]. Then n\ is fixed by the first eq. (10) and the values of the constants of motion take the 



values dictated by the choice of R and cr in (16) and (17), which change order by order in the 
expansion. 

Alternatively, one can use the parameters to improve the approximation of a particular 
orbit at fixed order in a. For example, one can require the constants of motion e and £ to 
have a fixed value at all orders in a from the first order onwards. However, in general this 
can not be done by adjusting the single parameter A„ at order n, and as a result the actual 
values of all parameters {R, o", Ai, A„) will have to be readjusted order by order. 

In practice, we prefer to use the freedom in the orbital parametrization to adjust all 
parameters so as to keep a growing number of orbital characteristics fixed order by order in 
the expansion: the orbital period at zeroth order, the radial co-ordinate of the periastron and 
the periastron advance at first order, the apastron co-ordinate at second order, etc., and to 
keep them fixed henceforth order by order in the expansion. Again, to this end the actual 
values of the parameters (i?, a, Ai, A„), and correspondingly the constants of motion e 
and i, have to be readjusted at every order n. For the second-order approximation presented 
above, these conditions read explicitl}{^ 

rpa = al + a[ + a2 + raa = aQ - a[ + + (18) 

whilst the total proper-time period Ar, observer time period At and angle Aip between 
periastra are 

ujAt = {ujo + (^^1 + = 2tt, (19) 

and hence 

ojAt = 27rao, uAip = 27ra^. (20) 
These four conditions determine the values of {R,a, Ai, A2) at second order. 

3 Numerical results 

The accuracy of our perturbation theory can be tested by comparison with the results of a 
purely numerical approach. As shown in the appendix of ref. p], we can parametrize any 
bound orbit in terms of a quasi semi- major axis a and a quasi eccentricity e such that 

l-'^O + ecosy). (21) 



l + ecosy{(f) \dtf J a 
The parameters (a, e) are determined by the constants of motion (e, i) via 



(22) 



of 
M 



+ r (3 + e^) = 0. 



^Our initial conditions are such that the test mass is at periastron at r = 0; therefore and a are 
necessarily negative. 
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In this section we perform a detailed comparison of both the first- and second-order pertur- 
bation expressions with a series of high-precision numerical solutions of the geodesic equation 
characterized by specific values of a and e. It is shown that the accuracy varies rather little 
with decreasing values of a, affirming that the epicycle approximations remain very accurate 
all the way up to the ISCO. 

As a typical example we consider an orbit with a = lOM and e = 0.1. For this orbit 
the radial co-ordinate of the periastron is rpa = 9.09091M, whilst the apastron is found at 
raa = ll.llllM. Recall, that the ISCO is the circular orbit at r = 6M. Furthermore, the 
angular shift between consecutive periastra is large: 5(p = Aip — 27r = 3.6561. This shows that 
for the orbit considered the effects of the curvature of space-time are large. In the following 
for convenience of numerical calculations we take M = 10 in some arbitrary units. In these 
units the proper time between the periastra is Ar = 2661.05, and the energy and angular 
momentum per unit of mass are 

e = 0.956568, i = 37.8235. (23) 

Now we construct the first-order epicycle approximation to this geodesic. From eqs. ([4|)-(|13|) 

we get 

rpa = R + a (Ai +nl) = R + ct^I- 



2tt „ 1 



At = — = 27rR . 



3M 
R 



M 1 - 6M ' (24) 



2tt ( 2,(7 Ai 1 



Ay = a^Ar = , ~" ( i_ ^ . 

These three conditions completely determine R, Ai and a, 

R = 101.274, Ai = 0.05047, a = -11.5699. (25) 

From these results we can compute the other dependent parameters and constants of motion: 

n[ = 0.845359, = 0.00236, e = eQ + aei = 0.956418, £ = 4 + ct^i = 37.8710. (26) 

Observe, that the first-order values of e and d. are accurate to about one per mille. We thus 
find the following explicit first-order approximation to the orbit: 

t(r) = 1.19347r + 11.998 sin(0.00236r), 

r(T) = 100.690- 9.78072 cos(0.00236t), (27) 

(^(r) = 0.00373 r + 0.3025 sin(0.00236 r) . 
Next we turn to the second order approximation to the exact orbit. As boundary conditions 
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Figure 1: The ratio ?'(T)/r„„„j(r) in first-order approximation (solid curve) and in second-order 
approximation (dashed curve). 




Figure 2: The difference 'p{t) — ipnum{T) in radians, in first-order approximation (solid curve) and in 
second-order approximation (dashed curve). 



we impose 



rpa = R + a (Ai + n[) + — (A2 + + m^) , 



raa = R + (T (Ai - n\) + — {A2-nl + , 



27r 2tt ( auj\ 

Ar = ~ — I 1 



(28) 



A93 



al At 



2-K 

UJqR 



R 



+ a\Rvl 



i^oRJ 



a'R 



2uji 



Given the values of periastron and apastron, and the periastron shift in angle and proper 
time as above, these equations determine R, Ai, A2 and cr; however, as the equations are 
quadratic in o", there are in general multiple solutions. We always choose the one with the 
smallest ratio of second-order to first-order contributions to <^(t); in practice we find that 
this also implies fastest convergence for r(r). Using this criterion our second-order solution 
is characterized by 



R = 100.046, Ai = -0.12613, A2 = 0.000228, a = -9.71137. 



(29) 



Observe, that the values of R, Ai and a differ from those at first order, as we are optimizing 
the fit of our approximation to a given orbit order-by-order in perturbation theory. For the 
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e = 0.1 


First order epicycles 


Second order epicycles 


a 


max. rel. difF. in r 


max. abs. diff. in tp 


max. rel. diff. in r 


max. abs. diff. in ip 


200 


1.05% 


0.018 


0.06% 


0.0015 


150 


1.03% 


0.020 


0.06% 


0.0018 


100 


1.08% 


0.029 


0.05% 


0.0040 


90 


1.09% 


0.024 


0.05% 


0.0070 


85 


1.04% 


0.042 


0.04% 


0.012 


80 


1.21% 


0.058 


0.04% 


0.026 


75 


1.47% 


0.028 


0.08% 


0.053 


70 


1.52% 


0.040 


0.12% 


0.012 


66 


2.20% 


0.022 


0.24% 


0.020 



Table 1: Accuracy of epicycle approximation for orbits of eccentricity e = 0.1 and various values of a, 
as listed in the first column. The second column presents the maximum relative differences between 
r(r) at first order and its numerical counterpart; the third column presents the maximum absolute 
difference between the azimuth fir) at first order and its numerical counterpart. The fourth and fifth 
column present the same information for the second-order approximations. 



other parameters the results (30) imply 

n[ = 1.02061, LJ = 0.00236, e = 0.95667, i = 37.8263. (30) 
Thus the explicit second-order solution reads 

t{T) = 1.19246T + 12.3903 sin(0.00236r) + 0.613326 sin(0.00472T), 



r(r) = 101.378 - 10.101 cos(0.00236 r) - 0.368364 cos(0.00472 r), (31) 

p{t) = 0.0037351 r + 0.312458 sin(0.00236 r) + 0.0174545 sin(0.00472 r) . 

To compare these results with the purely numerical solution in figs. 1 and 2 we plot the ratio of 
our perturbative ?'('?") to the numerical result fnum 

(r), and the difference of our perturbative 
(P{t) and the numerical values (pnum{T). These figures show, that our approximations to the 
radial distance are accurate to better than 1.1 % at first order and to better than 0.05 % at 
second order, whilst the accuracy of the azimuthal co-ordinate is 0.029 radians at first order, 
and 0.004 radians at second order. 

To test the accuracy of our perturbation theory in the regime of strong curvature we have 
performed a series of similar calculations for values of a ranging from 20M to 6.6M, all with 
quasi eccentricity e = 0.1, such that the smallest orbit actually reaches the ISCO at 6M at 
its periastron. Again, for numerical purposes we have chosen to set the central mass equal to 
M = 10 in arbitrary units. The comparison with numerical calculations of the corresponding 
orbits is summarized in table 1. 
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4 Gravitational waves 



Having in hand the perturbative solution of the geodesic equations for test masses in a 
Schwarzschild background, we can use them to compute the gravitational radiation emit- 
ted by extreme mass ratio binary systems in which the larger mass is non-rotating, at least 
in the limit in which radiation reaction effects can be ignored. The two linearized fluctuation 
modes of the Schwarzschild geometry are described by the Regge- Wheeler and Zerilli-Moncrief 
equations |im [TT\ \n\ . A test mass moving in a Schwarzschild background acts as a source for 
such fluctuations, giving rise to a specific type of source terms in the fluctuation equations. 

In a multipole expansion the fluctuation equations for the amplitudes X = {RW, ZM) 
denoting the Regge- Wheeler or Zerilli-Moncrief modes, take the form 

{Us - Vk{r)) ¥^^\t,r) = Sf{t,r), (32) 

where the radial d'Alembert operator is given by 

□5 = -df + fdrfdr, fir) = 1 - — . (33) 

r 

The potentials V^(r) depend on the radial co-ordinate r and the mode-index I: 



fu, 6Af 36M2 72M3\ 



(34) 



where A = (/ + 2) (/ — 1). Using the results of section [2] the source terms can be computed 
from the energy-momentum tensor for a point mass ^ 

Tf^u = j -^=Ufj,Uy5^ {x - Xp{T))dT, (35) 

to take the form 

5^ = Gf{r,{t)) Sir - r,{t)) + Fjr(r, r,{t)) 5'{r - T,{t)), (36) 

in which 

^RW\^^^vV^)) - ^(/ + 2)! r 'ff ' 

Gnw{r,r,[t)) - 32vr^^^^/(rp)|^l-— j — H^^^ +^^^ — — , 

(37) 



and 



(A + 2)(A+^j ^ 
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2/(rp) 



(A + 2) A + 



rp J 



(A + 2) A + ^^ 



6M'\ U* 



(A + 2) A + 



6M 
rp 



+ 



U U 



647r/i /(rp) 
/(/ + 1) _^ 6JU^ 



*lrn 



(/-2)! 



(38) 



Here y'"^, (X^,Z^) and {U^ab^V^.W^^) represent standard scalar, vector and tensor 
harmonics [161 I15j . whilst the four- velocity components are to be taken from sect. [2j 
Details will be given in |21j . 



The fluctuation equations (32) have been worked out and solved in a fully numerical 

We have numerically solved 



(39) 



approach for the source terms and fluctuations in refs. [10j-[T5l 
the equations for the asymptotic gravitational- wave amplitude at distance r: 



h+{t) - ihy,{t) 



starting from the fluctuation equations in full analytic form, using the algorithm of Lousto 



and Price [13], for the second-order approximation to the orbit discussed above, eqs. (29)-(31 ). 



From the amplitudes one can calculate the power emitted in terms of gravitational waves 
from the Regge- Wheeler and Zerilli-Moncrief functions by evaluating the expression |14^ [T5] 



P 



1 



(Z + 2)! 



647r ^ (I 

Im 



2)\ 



Im 
ZM 



+ 4 



1- 



Im 
RW 



(40) 



where the overdot denotes a derivative w.r.t. co-ordinate time t. Also, the angular momentum 
per unit time emitted by gravitational waves follows as 



dL 

~dt 



1287r 



E 

Im 



m{l + 2)\ 



^*Ryv{t')dt'] +C.C. 



(41) 



in which c.c. stand for the complex conjugate. 

We have performed the computation of the gravitational-wave power for the full series of 
orbits listed in table [T| which presents the power averaged over an orbital revolution. The 
results in units of the mass ratio {M/ fi)"^ are presented in the second column of table [2j 
They can be compared to the power as calculated by the Peters-Mathews equation in the 
Newtonian approximation |20j in the third column. To reach the required accuracy it is 
sufficient to compute the modes up to / = 4, m = 4 and add them as indicated in eqs. (39) 



- (41); higher-order contributions do not change the results. In fact, the main contribution 
comes from the modes with Z = 2, as is to be expected from the quadrupole nature of free 
gravitational waves. The numerical accuracy has furthermore been checked by requiring the 
amplitude and power to remain the same upon reducing the grid scale by a factor two. 

As an example, the I = 2 modes are shown for the reference orbit with a = lOM and 
e = 0.1 in figs. [3] and m The average power emitted in this orbit in the second-order epicycle 
approximation is 



(P) = 6.294 X 10 



-5 



(42) 
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a 


{-^ / epicycle 


{P)PM 


1 1 • rr 

rei. diii. 


200 


2.033 10"*^ 


2.033 lO""^ 


+0.2% 


150 


8.139 10"*^ 


8.555 lO"'* 


-5.1% 


100 


6.294 10"^ 


6.496 10"^ 


-3.2% 


90 


1.083 10-^ 


1.100 10-4 


-1.6% 


85 


1.475 10-^ 


1.464 10-4 


-0.5% 


80 


1.991 10-^ 


1.983 10-4 


+0.4% 


75 


2.944 10-^ 


2.738 10-4 


+7.6% 


70 


4.209 10"^ 


3.865 10-4 


+8.9% 


66 


5.869 10-4 


5.187 10-4 


+13% 



Table 2: Average power (P) emitted in gravitational waves by a system of a star of mass in a bound 
orbit of eccentricity e = 0.1 around a black hole of mass M = 10 for various values of a, as listed in 
the first column. The second column presents the power in the second order epicycle approximation, 
whereas the third column presents the power as calculated by the Peters-Mathews equation. The third 
column presents the relative difference in percents. All powers are computed to a numerical accuracy 
of 0.1% and are stated in units A'P / fi^ . 



which differs a mere 0.3 % from the fully numerical result |18^ I19j 

(P) = 6.318 X 10"^ (43) 




Figure 3: Real part of the Z = 2, m = 2 Zerilli-Moncrief function (left) and power emitted by this 
mode (right) as a function of time t, for an observer located at distance r — 500A/. 

Finally, we proceed to calculate the emitted power and angular momentum per unit time 
for a series of eccentric orbits that are adiabatically related to each other by the emission 
of gravitational radiation. That is, we calculate the average power and angular momentum 
emitted by gravitational waves for each given eccentric orbit, and use these numbers to update 
the values of e and ^ in a discrete step; the newly found values for e and £ then correspond to 
the next orbit in the series. The discrete step is chosen as follows: the next orbit will always 
be chosen to be the one that has a value i that is 1% smaller than that of the current orbit. 
In this way, we find in practice that successive orbits have a percentual change in periastra 
of less than 4%, which justifies the adiabatic approximation. The results are shown in Table 
[sj As before, the size of the grid was chosen such that the values for < P > and < L > do 
not change more than at the 0.1% level when taking a grid size a factor of two smaller still. 
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Figure 4: Real part of the I = 2, m — 1 Regge- Wheeler function (left) and power emitted by this 
mode (right) as a function of time t, for an observer located at distance r = 500A/. 

As can be seen, under influence of the emission of gravitational waves the eccentricity 
of the orbits decreases, and the orbits become more circular after every discrete step. As a 
result the epicycle expansion is expected to become increasingly accurate for successive orbits, 
and indeed this is seen to be the case: the radial orbital function r becomes almost twice 
as accurate in the course of the orbits considered. In contrast, the precision of the angular 
coordinate does not improve, but remains more or less fixed around the value of 0.008 
radians. This conforms to expectations, as in the second order epicycle expansion we use two 
of the four boundary conditions to fix he time and angular shift between successive periastra, 
leaving only two boundary conditions to fit the orbital functions r and (p. In the calculation 
presented, these remaining two boundary conditions were used to fix the values of the radial 
positions of the periastron and apastron, which in practice also makes the orbital function if 
very accurate, but without causing the accuracy to increase with decreasing eccentricity. 



5 Discussion and conclusions 

Comparing the results of the epicycle approximation with the Peters-Mathews calculations as 
collected in table [2j one observes that both at large distances and close to the ISCO the power 
computed by our relativistic procedure exceeds that of the newtonian approximation, whilst 
in contrast for intermediate distances the relativistic power is lower. This can be understood 
by the opposite effect of two factors: on the one hand the precession of the periastron shows 
that the orbital velocity and acceleration in the relativistic orbit is higher than that in the 
corresponding Kepler orbit; on the other hand in the relativistic computation the redshift of 
the gravitational waves lowers the power emitted as measured by a distant observer. 

Comparison with purely numerical calculations shows that for the series of orbits consid- 
ered the second-order epicycle approximation is very accurate, at the level of one part in a 
thousand. The agreement is less for orbits with larger eccentricity; for e = 0.2 the agreement 
is still very good for large orbits such as a = 20M, but closer to the ISCO the accuracy 
becomes of the order of 1-5%. To improve on this it is necessary to include the third-order 
epicycle contribution. Indeed, the source of deviations is the dependence on the azimuth an- 
gle ^p{t), which appears as argument in the tensor spherical harmonics. The angular velocity 
can be improved significantly without compromizing the radial accuracy only by including 
the third-order epicycle terms [1]. However, as the emission of gravitational radiation tends 
to lead to significant loss of angular momentum, thus decreasing the eccentricity in the last 
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a 


e 


<P> 


< L> 


Aa/a 


Ae/e 


At 


max.rel.diff. r 


max.abs.difF. ip 


120.0 


0.1500 


2.536 


0.9961 


n.a. 


n.a. 


n.a. 


0.20% 


0.0085 


116.4 


0.1447 


2.951 


1.125 


3.1% 


3.7% 


4.021 


0.18% 


0.0080 


112.8 


0.1410 


3.451 


1.258 


3.2% 


2.6% 


3.526 


0.16% 


0.0080 


109.3 


0.1378 


4.051 


1.410 


3.2% 


2.4% 


3.121 


0.15% 


0.0081 


105.8 


0.1349 


4.785 


1.589 


3.3% 


2.2% 


2.756 


0.14% 


0.0083 


102.3 


0.1322 


5.683 


1.794 


3.4% 


2.0% 


2.423 


0.13% 


0.0087 


98.7 


0.1300 


6.774 


2.039 


3.6% 


1.8% 


2.123 


0.12% 


0.0094 


95.2 


0.1287 


8.193 


2.326 


3.7% 


1.0% 


1.850 


0.11% 


0.011 


91.6 


0.1273 


9.853 


2.661 


3.9% 


1.1% 


1.604 


0.11% 


0.013 



Table 3: A series of eccentric orbits parametrized by a and e, as listed in the first two columns. 
The orbits are related by the emission of gravitational waves in an adiabatic way, as explained in the 
main text. The average power emitted is given in the third column in imits 10^^"'' (M//i)^, and the 
average angular momentum emitted per unit time is given in the fourth column in units 10^^^ (M//x^] . 
The fifth and sixth columns list the percentual change in a and e compared to the next larger orbit, 
showing explicitly the inspiral and circularization due to the emission of gravitational waves. The 
seventh column presents the Schwarzschild time taken to make the discrete step from the previous 
orbit to the current, in units 10^^ seconds. The eighth column presents the maximum relative difference 
between the second-order epicycle radial orbital function r(T) and its numerical counterpart; the last 
column presents the maximum absolute difference between the second-order epicycle angular orbital 
function (^(r) and its numerical counterpart. The mass of the black hole and that of the companion 
star are, in arbitrary units, set to M = 10 and /x = 10~*. 



stage of inspiral, in practice the second-order epicycle approximation leads mostly to very 
acceptable results. 
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